JSDM models test
Data description
Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
me5 <- load_object("model-2019-04-09-19-02-16.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me5)## Summary for model '/var/folders/2z/0sxx_8ts2pxcp028sy41wfgc0000gn/T//RtmpeNqPAL/file1219c690e12c'
## Saved parameters: B Rho EnvRho Tau
## MCMC ran in parallel for 34.541 minutes at time 2019-04-09 18:27:42.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.#me5$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me5$R$B),4)))## Maximum Rhat value for Beta: 1.136cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me5$R$Rho),4)))## Maximum Rhat value for Rho: NAcat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me5$R$EnvRho),4)))## Maximum Rhat value for EnvRho: 1.1004cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me5$R$Tau),4)))## Maximum Rhat value for Tau: 1.4694me5$n.eff## $B
## [,1] [,2] [,3]
## [1,] 29 155 32
## [2,] 257 547 283
## [3,] 756 1000 1000
## [4,] 103 593 94
## [5,] 181 204 177
##
## $Rho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 86 66 284 232
## [2,] 86 1 132 204 350
## [3,] 66 132 1 195 53
## [4,] 284 204 195 1 195
## [5,] 232 350 53 195 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 54 39 53 48
## [2,] 54 1 231 129 277
## [3,] 39 231 1 189 273
## [4,] 53 129 189 1 181
## [5,] 48 277 273 181 1
##
## $Tau
## [,1] [,2] [,3] [,4] [,5]
## [1,] 12 167 129 116 376
## [2,] 167 29 411 82 216
## [3,] 129 411 33 491 159
## [4,] 116 82 491 21 251
## [5,] 376 216 159 251 59me5$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000Environment filtering 5 species
JSDM
Model for \(\textbf{EnvEvenSp5}\)
data<-sim_data$EnvEvenSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
Y_cor<-cor(data$Y)
to_prec<-function(m){
n<-dim(m)[1]
Tau_n<-matrix(nrow=n, ncol=n)
for (j in 1:n) {
for (k in 1:n){
Tau_n[j, k] <- -m[j, k]/sqrt((m[j,j]*m[k,k]))
}
}
return(Tau_n)
}
#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Tau signif")
GJAM
##to setup chains parameters
data<-sim_data$EnvEvenSp5
#fit_gjam(data,2000,1000,"./gjam_models/gjam5env.rda")
load_gjam(data,name="./gjam_models/gjam5env.rda")##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 9.5 3.13 5.08 16.4
## env2 31.9 7.93 20.50 50.6
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05
## intercept -0.0312 0.646 1.2500 0.704 -1.200
## env -2.4200 -1.800 0.0438 1.710 2.870
## env2 0.9440 -0.527 -0.9410 -0.598 -0.665
## RMSPE 0.2780 0.327 0.3760 0.309 0.327
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.0312 0.0775 -0.1940 0.107
## sp01_env -2.4200 0.1170 -2.6500 -2.200 *
## sp01_env2 0.9440 0.0938 0.7780 1.130 *
## sp02_intercept 0.6460 0.0651 0.5220 0.771 *
## sp02_env -1.8000 0.0911 -1.9900 -1.630 *
## sp02_env2 -0.5270 0.0741 -0.6700 -0.389 *
## sp03_intercept 1.2500 0.0672 1.1100 1.380 *
## sp03_env 0.0438 0.0465 -0.0447 0.132
## sp03_env2 -0.9410 0.0764 -1.0900 -0.781 *
## sp04_intercept 0.7040 0.0865 0.5590 0.881 *
## sp04_env 1.7100 0.0897 1.5500 1.890 *
## sp04_env2 -0.5980 0.0608 -0.7170 -0.486 *
## sp05_intercept -1.2000 0.0953 -1.3700 -1.000 *
## sp05_env 2.8700 0.1380 2.6300 3.170 *
## sp05_env2 -0.6650 0.0610 -0.7840 -0.545 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 5 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.325, and the DIC is 25294. Computation involved 2000 Gibbs steps, with a burnin of 1000.
HMSC
#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="F_t",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
if (label=="F_t"){
Y_data = subset(data, select = -env)
ns<- ncol(Y_data)
np <- nrow(Y_data)
X<-scale(poly(data$env[1:np], 2))
colnames(X)<-c("env","env2")
studyDesign = data.frame(sample = as.factor(1:np))
rL = HmscRandomLevel(units = studyDesign$sample)
m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
studyDesign = studyDesign, ranLevels = list(sample = rL))
m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
save(m, file = name)
return(m)
}
if (label=="L_d"){
return(load_object(name))
}
}
data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")Convergence:
hm_conv<-function(mod){
codaList = convertToCodaObject(mod)
#convergence histograms
hist(effectiveSize(codaList$Beta), main="ess(beta)")
hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf, main="psrf(beta)")
hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)")
hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf, main="psrf(omega)")
}
hm_conv(hm_mod)
Study of interactions
hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
getOmega = function(a,r=1)
return(crossprod(a$Lambda[[r]]))
ns<-mod$ns
postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
postOmega<-abind(postOmega1,postOmega2,along=3)
postOmegaMean = apply(postOmega,c(1,2),mean)
postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
postR<-array(dim=c(ns,ns,nchains*nsamples))
for(i in 1:dim(postOmega)[3])
postR[,,i]<-stats::cov2cor(postOmega[,,i])
postRMean = apply(postR,c(1,2),mean)
postRUp=apply(postR,c(1,2),quantile,0.95)
postRLo=apply(postR,c(1,2),quantile,0.05)
Tau = solve(postOmegaMean)
Tau_n = cov2cor(Tau)
Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
# Omegacor<- computeAssociations(m)
# supportLevel<- 0.95
# toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
# corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(hm_mod$Y), diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Correlation cor(Y)")
corrplot(postRMean, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("R")
corrplot(Toplot_R, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Plot only non zero value")
corrplot(Tau_n, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Partial correlation matrix")
corrplot(interact, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("True interactions")
}
hm_inter(hm_mod)
Environmental filtering 10 species
JSDM
me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.#me10$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me10$R$B),4)))## Maximum Rhat value for Beta: 1.5602cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me10$R$Rho),4)))## Maximum Rhat value for Rho: NAcat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me10$R$EnvRho),4)))## Maximum Rhat value for EnvRho: 1.7953cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me10$R$Tau),4)))## Warning in max(me10$R$Tau): no non-missing arguments to max; returning -Inf## Maximum Rhat value for Tau: -Infme10$n.eff## $B
## [,1] [,2] [,3]
## [1,] 34 51 38
## [2,] 33 121 31
## [3,] 130 137 82
## [4,] 271 1000 645
## [5,] 824 593 1000
## [6,] 756 670 1000
## [7,] 729 670 259
## [8,] 58 69 61
## [9,] 45 47 47
## [10,] 11 47 10
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 119 49 54 84 102 66 519 121 186
## [2,] 119 1 162 121 81 136 184 283 224 219
## [3,] 49 162 1 755 53 34 120 72 256 217
## [4,] 54 121 755 1 357 80 82 75 22 273
## [5,] 84 81 53 357 1 269 76 84 62 59
## [6,] 102 136 34 80 269 1 249 404 136 78
## [7,] 66 184 120 82 76 249 1 391 363 30
## [8,] 519 283 72 75 84 404 391 1 36 66
## [9,] 121 224 256 22 62 136 363 36 1 157
## [10,] 186 219 217 273 59 78 30 66 157 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 60 74 80 59 39 39 34 89 18
## [2,] 60 1 36 56 44 49 39 86 230 9
## [3,] 74 36 1 144 134 112 110 51 249 23
## [4,] 80 56 144 1 505 379 214 67 74 14
## [5,] 59 44 134 505 1 632 1000 113 106 16
## [6,] 39 49 112 379 632 1 801 103 86 14
## [7,] 39 39 110 214 1000 801 1 128 102 14
## [8,] 34 86 51 67 113 103 128 1 42 17
## [9,] 89 230 249 74 106 86 102 42 1 23
## [10,] 18 9 23 14 16 14 14 17 23 1me10$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000data<-sim_data$EnvEvenSp10
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Correlation cor(Y)")
corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho")
corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho signif")
corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho")
corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("True interactions")
#corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
#title("Tau")
#corrplot(Tau_n*(!me10$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
#title("Tau")
}
plot_cor_jsdm(me10,data$Y)
Gjam
data<-sim_data$EnvEvenSp10
#fit_gjam(data,5000,500,"./gjam_models/gjam10env.rda")
load_gjam(data,name="./gjam_models/gjam10env.rda")##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 233.0 58.0 120 354.0
## env2 56.7 11.8 37 82.9
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08 sp09
## intercept -0.2060 0.3800 1.310 0.889 1.020 0.742 0.638 0.140 0.332
## env -2.0500 -2.5700 -2.300 -1.030 -0.189 0.542 1.090 1.740 2.520
## env2 0.0656 -0.0531 0.338 -0.533 -0.805 -0.838 -0.690 -0.425 0.640
## RMSPE 0.3300 0.2980 0.317 0.371 0.414 0.421 0.374 0.345 0.298
## sp10
## intercept -0.577
## env 2.370
## env2 0.677
## RMSPE 0.276
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.2060 0.0670 -0.33700 -0.0733 *
## sp01_env -2.0500 0.1190 -2.29000 -1.8300 *
## sp01_env2 0.0656 0.0992 -0.14100 0.2380
## sp02_intercept 0.3800 0.1030 0.18000 0.5510 *
## sp02_env -2.5700 0.1370 -2.83000 -2.2900 *
## sp02_env2 -0.0531 0.0790 -0.19700 0.1120
## sp03_intercept 1.3100 0.0829 1.15000 1.4700 *
## sp03_env -2.3000 0.1260 -2.55000 -2.0600 *
## sp03_env2 0.3380 0.0712 0.20500 0.4860 *
## sp04_intercept 0.8890 0.0585 0.77400 1.0000 *
## sp04_env -1.0300 0.0659 -1.16000 -0.9070 *
## sp04_env2 -0.5330 0.0580 -0.64200 -0.4160 *
## sp05_intercept 1.0200 0.0621 0.90800 1.1500 *
## sp05_env -0.1890 0.0488 -0.28500 -0.0949 *
## sp05_env2 -0.8050 0.0759 -0.96000 -0.6590 *
## sp06_intercept 0.7420 0.0635 0.61600 0.8650 *
## sp06_env 0.5420 0.0558 0.43200 0.6460 *
## sp06_env2 -0.8380 0.0634 -0.96100 -0.7130 *
## sp07_intercept 0.6380 0.0795 0.50600 0.8160 *
## sp07_env 1.0900 0.0805 0.93900 1.2500 *
## sp07_env2 -0.6900 0.0595 -0.80300 -0.5700 *
## sp08_intercept 0.1400 0.0759 0.00798 0.2950 *
## sp08_env 1.7400 0.1240 1.48000 1.9600 *
## sp08_env2 -0.4250 0.1000 -0.65400 -0.2750 *
## sp09_intercept 0.3320 0.0539 0.22600 0.4370 *
## sp09_env 2.5200 0.1660 2.22000 2.8500 *
## sp09_env2 0.6400 0.0700 0.50400 0.7780 *
## sp10_intercept -0.5770 0.0696 -0.70200 -0.4380 *
## sp10_env 2.3700 0.1130 2.16000 2.6000 *
## sp10_env2 0.6770 0.0820 0.49900 0.8200 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 10 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.348, and the DIC is 114934. Computation involved 5000 Gibbs steps, with a burnin of 500.
#gje10<-load_object("./gjam_models/gjam10env.rda")
#to check posterior density of s in Sigma
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))HMSC
data<-sim_data$EnvEvenSp10
#hm_mod<-fit_hmsc(data,"F_t",name="./HMmodels/hm10env.rda" )
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm10env.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod)
Environmental filtering 20 species
JSDM
me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.#me20$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me20$R$B),4)))## Maximum Rhat value for Beta: 1.2174cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me20$R$Rho),4)))## Maximum Rhat value for Rho: NAcat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me20$R$EnvRho),4)))## Maximum Rhat value for EnvRho: 1.3236cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me20$R$Tau),4)))## Warning in max(me20$R$Tau): no non-missing arguments to max; returning -Inf## Maximum Rhat value for Tau: -Infme20$n.eff## $B
## [,1] [,2] [,3]
## [1,] 114 142 148
## [2,] 39 89 30
## [3,] 21 99 21
## [4,] 23 30 22
## [5,] 26 46 30
## [6,] 88 141 85
## [7,] 39 35 42
## [8,] 675 220 85
## [9,] 1000 265 1000
## [10,] 674 700 218
## [11,] 598 276 1000
## [12,] 180 413 334
## [13,] 438 378 386
## [14,] 366 129 155
## [15,] 56 44 83
## [16,] 53 97 63
## [17,] 23 45 26
## [18,] 42 59 30
## [19,] 19 68 21
## [20,] 37 116 30
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 181 118 141 73 46 71 89 138 25 49 63 247
## [2,] 181 1 78 200 111 18 147 186 58 40 121 74 37
## [3,] 118 78 1 76 43 28 20 59 124 212 85 91 74
## [4,] 141 200 76 1 354 143 269 122 173 42 66 90 29
## [5,] 73 111 43 354 1 533 128 100 189 370 176 57 42
## [6,] 46 18 28 143 533 1 661 73 175 97 963 78 68
## [7,] 71 147 20 269 128 661 1 212 143 324 100 84 125
## [8,] 89 186 59 122 100 73 212 1 107 395 156 115 115
## [9,] 138 58 124 173 189 175 143 107 1 813 131 143 109
## [10,] 25 40 212 42 370 97 324 395 813 1 106 67 49
## [11,] 49 121 85 66 176 963 100 156 131 106 1 431 244
## [12,] 63 74 91 90 57 78 84 115 143 67 431 1 145
## [13,] 247 37 74 29 42 68 125 115 109 49 244 145 1
## [14,] 567 119 24 413 41 348 26 353 197 135 121 178 318
## [15,] 111 174 59 353 54 31 26 295 42 48 1000 66 75
## [16,] 77 136 59 56 76 47 402 80 44 101 245 450 53
## [17,] 105 82 67 58 70 171 44 80 51 43 57 88 30
## [18,] 316 422 58 81 425 152 110 131 55 164 143 72 30
## [19,] 63 212 130 112 125 110 199 78 56 29 81 45 31
## [20,] 143 364 127 489 821 208 563 80 273 74 318 679 58
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 567 111 77 105 316 63 143
## [2,] 119 174 136 82 422 212 364
## [3,] 24 59 59 67 58 130 127
## [4,] 413 353 56 58 81 112 489
## [5,] 41 54 76 70 425 125 821
## [6,] 348 31 47 171 152 110 208
## [7,] 26 26 402 44 110 199 563
## [8,] 353 295 80 80 131 78 80
## [9,] 197 42 44 51 55 56 273
## [10,] 135 48 101 43 164 29 74
## [11,] 121 1000 245 57 143 81 318
## [12,] 178 66 450 88 72 45 679
## [13,] 318 75 53 30 30 31 58
## [14,] 1 111 183 58 75 52 86
## [15,] 111 1 221 274 60 333 59
## [16,] 183 221 1 128 70 469 101
## [17,] 58 274 128 1 233 118 159
## [18,] 75 60 70 233 1 117 161
## [19,] 52 333 469 118 117 1 172
## [20,] 86 59 101 159 161 172 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 99 33 72 38 483 109 72 139 273 137 356 205
## [2,] 99 1 68 111 59 42 33 28 34 36 43 44 50
## [3,] 33 68 1 51 48 33 47 31 30 32 41 28 39
## [4,] 72 111 51 1 36 51 22 27 40 40 60 49 58
## [5,] 38 59 48 36 1 32 55 57 50 35 66 42 69
## [6,] 483 42 33 51 32 1 58 70 102 233 76 227 133
## [7,] 109 33 47 22 55 58 1 65 61 50 67 62 74
## [8,] 72 28 31 27 57 70 65 1 228 224 195 178 357
## [9,] 139 34 30 40 50 102 61 228 1 351 520 689 1000
## [10,] 273 36 32 40 35 233 50 224 351 1 169 794 489
## [11,] 137 43 41 60 66 76 67 195 520 169 1 310 1000
## [12,] 356 44 28 49 42 227 62 178 689 794 310 1 856
## [13,] 205 50 39 58 69 133 74 357 1000 489 1000 856 1
## [14,] 312 74 32 37 48 149 91 167 351 214 287 456 201
## [15,] 62 38 43 94 58 90 68 152 184 131 103 114 112
## [16,] 148 62 27 315 28 174 39 56 128 121 76 211 132
## [17,] 58 59 56 30 32 90 20 29 36 56 29 47 48
## [18,] 44 34 25 53 64 42 34 52 77 39 123 57 108
## [19,] 100 37 127 42 53 33 50 23 25 29 30 32 31
## [20,] 46 18 15 29 20 89 35 34 56 47 45 73 68
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 312 62 148 58 44 100 46
## [2,] 74 38 62 59 34 37 18
## [3,] 32 43 27 56 25 127 15
## [4,] 37 94 315 30 53 42 29
## [5,] 48 58 28 32 64 53 20
## [6,] 149 90 174 90 42 33 89
## [7,] 91 68 39 20 34 50 35
## [8,] 167 152 56 29 52 23 34
## [9,] 351 184 128 36 77 25 56
## [10,] 214 131 121 56 39 29 47
## [11,] 287 103 76 29 123 30 45
## [12,] 456 114 211 47 57 32 73
## [13,] 201 112 132 48 108 31 68
## [14,] 1 80 101 40 47 33 48
## [15,] 80 1 1000 63 127 19 172
## [16,] 101 1000 1 96 337 27 188
## [17,] 40 63 96 1 33 355 35
## [18,] 47 127 337 33 1 17 42
## [19,] 33 19 27 355 17 1 18
## [20,] 48 172 188 35 42 18 1me20$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000data<-sim_data$EnvEvenSp20
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm(me20,data$Y)
Gjam
data<-sim_data$EnvEvenSp20
#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
load_gjam(data,name="./gjam_models/gjam20env.rda")##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 335 61.7 222 461
## env2 157 28.3 110 219
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08 sp09
## intercept -0.533 0.0688 0.500 0.144 0.327 0.498 0.720 0.896 0.906
## env -1.850 -2.5200 -2.480 -2.020 -1.820 -1.760 -1.220 -0.753 -0.516
## env2 0.187 0.3400 0.451 -0.279 -0.507 -0.493 -0.651 -0.641 -0.815
## RMSPE 0.337 0.3000 0.309 0.331 0.332 0.333 0.365 0.398 0.406
## sp10 sp11 sp12 sp13 sp14 sp15 sp16 sp17 sp18
## intercept 1.130 1.120 1.050 0.689 0.610 0.435 0.440 0.113 -0.07100
## env -0.200 0.133 0.551 0.870 1.140 1.730 2.460 1.990 2.35000
## env2 -0.931 -0.855 -0.860 -0.794 -0.671 -0.452 -0.227 -0.271 0.00271
## RMSPE 0.398 0.411 0.386 0.399 0.363 0.339 0.308 0.330 0.30300
## sp19 sp20
## intercept -0.0697 -0.329
## env 2.5800 2.190
## env2 0.3120 0.530
## RMSPE 0.3030 0.318
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.53300 0.1340 -0.764000 -0.3250 *
## sp01_env -1.85000 0.1340 -2.100000 -1.5800 *
## sp01_env2 0.18700 0.0717 0.055600 0.3320 *
## sp02_intercept 0.06880 0.0972 -0.095900 0.2760
## sp02_env -2.52000 0.1560 -2.810000 -2.2200 *
## sp02_env2 0.34000 0.0720 0.206000 0.4840 *
## sp03_intercept 0.50000 0.1040 0.295000 0.6900 *
## sp03_env -2.48000 0.1360 -2.770000 -2.2400 *
## sp03_env2 0.45100 0.0813 0.310000 0.6170 *
## sp04_intercept 0.14400 0.0783 0.009940 0.3210 *
## sp04_env -2.02000 0.0913 -2.200000 -1.8400 *
## sp04_env2 -0.27900 0.0825 -0.432000 -0.1160 *
## sp05_intercept 0.32700 0.0843 0.172000 0.4850 *
## sp05_env -1.82000 0.1010 -2.020000 -1.6200 *
## sp05_env2 -0.50700 0.0716 -0.653000 -0.3720 *
## sp06_intercept 0.49800 0.0786 0.342000 0.6400 *
## sp06_env -1.76000 0.0786 -1.920000 -1.6100 *
## sp06_env2 -0.49300 0.0669 -0.631000 -0.3680 *
## sp07_intercept 0.72000 0.0628 0.598000 0.8430 *
## sp07_env -1.22000 0.0844 -1.370000 -1.0500 *
## sp07_env2 -0.65100 0.0816 -0.819000 -0.4980 *
## sp08_intercept 0.89600 0.0702 0.766000 1.0400 *
## sp08_env -0.75300 0.0678 -0.886000 -0.6230 *
## sp08_env2 -0.64100 0.0614 -0.759000 -0.5180 *
## sp09_intercept 0.90600 0.0712 0.776000 1.0500 *
## sp09_env -0.51600 0.0835 -0.663000 -0.3500 *
## sp09_env2 -0.81500 0.0669 -0.942000 -0.6790 *
## sp10_intercept 1.13000 0.0655 1.000000 1.2600 *
## sp10_env -0.20000 0.0573 -0.305000 -0.0862 *
## sp10_env2 -0.93100 0.0676 -1.060000 -0.8010 *
## sp11_intercept 1.12000 0.0684 0.987000 1.2500 *
## sp11_env 0.13300 0.0565 0.022800 0.2410 *
## sp11_env2 -0.85500 0.0630 -0.978000 -0.7310 *
## sp12_intercept 1.05000 0.0730 0.906000 1.2000 *
## sp12_env 0.55100 0.0747 0.420000 0.7010 *
## sp12_env2 -0.86000 0.0752 -1.010000 -0.7180 *
## sp13_intercept 0.68900 0.0591 0.570000 0.8000 *
## sp13_env 0.87000 0.0813 0.709000 1.0300 *
## sp13_env2 -0.79400 0.0668 -0.922000 -0.6560 *
## sp14_intercept 0.61000 0.0641 0.481000 0.7300 *
## sp14_env 1.14000 0.0675 1.020000 1.2800 *
## sp14_env2 -0.67100 0.0665 -0.800000 -0.5470 *
## sp15_intercept 0.43500 0.0588 0.325000 0.5530 *
## sp15_env 1.73000 0.1250 1.520000 1.9800 *
## sp15_env2 -0.45200 0.0641 -0.579000 -0.3340 *
## sp16_intercept 0.44000 0.0521 0.337000 0.5420 *
## sp16_env 2.46000 0.1050 2.260000 2.6700 *
## sp16_env2 -0.22700 0.0501 -0.323000 -0.1280 *
## sp17_intercept 0.11300 0.0563 0.000618 0.2230 *
## sp17_env 1.99000 0.1310 1.760000 2.2600 *
## sp17_env2 -0.27100 0.0645 -0.394000 -0.1480 *
## sp18_intercept -0.07100 0.0549 -0.181000 0.0357
## sp18_env 2.35000 0.1040 2.140000 2.5500 *
## sp18_env2 0.00271 0.0659 -0.130000 0.1260
## sp19_intercept -0.06970 0.0875 -0.241000 0.0747
## sp19_env 2.58000 0.1220 2.350000 2.8300 *
## sp19_env2 0.31200 0.0517 0.216000 0.4120 *
## sp20_intercept -0.32900 0.0556 -0.433000 -0.2160 *
## sp20_env 2.19000 0.0953 2.010000 2.3800 *
## sp20_env2 0.53000 0.0616 0.409000 0.6490 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 20 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.351, and the DIC is 371075. Computation involved 5000 Gibbs steps, with a burnin of 500.
#gje20<-load_object("./gjam_models/gjam20env.rda")
#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))Gjam dimension reduction
data<-sim_data$EnvEvenSp20
#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
#load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")
#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
rl <- list(r = 8, N = 20)
formula<-as.formula( ~env+ env2)
ml <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
####fit
mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)##
## Observations and responses:
## [1] 500 20## Warning in .setupReduct(modelList, S, Q, n): dimension reduction requires
## reductList$N < no. responses##
## Dimension reduced from 20 X 20 -> 20 X 8 responses
## ===========================================================================## expanding covariance chains## =========================================================================== save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
summary(mod_gjam1)##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 0.963 0.0606 0.844 1.08
## env2 1.310 0.0698 1.180 1.45
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08
## intercept -0.344 -0.1010 0.00677 0.111 0.292 0.342 0.469 0.483
## env -0.473 -0.4810 -0.48100 -0.475 -0.470 -0.463 -0.430 -0.356
## env2 0.276 0.0944 0.00263 -0.110 -0.262 -0.329 -0.441 -0.460
## RMSPE 0.406 0.4050 0.40800 0.410 0.403 0.402 0.408 0.424
## sp09 sp10 sp11 sp12 sp13 sp14 sp15 sp16 sp17
## intercept 0.484 0.488 0.488 0.487 0.478 0.475 0.392 0.276 0.177
## env -0.362 -0.091 0.198 0.462 0.423 0.457 0.466 0.476 0.478
## env2 -0.478 -0.484 -0.475 -0.482 -0.478 -0.454 -0.347 -0.259 -0.171
## RMSPE 0.431 0.441 0.446 0.427 0.425 0.405 0.406 0.407 0.408
## sp18 sp19 sp20
## intercept -0.000261 -0.0922 -0.277
## env 0.481000 0.4800 0.471
## env2 0.009230 0.1010 0.246
## RMSPE 0.404000 0.4060 0.403
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.344000 0.0705 -0.4720 -0.2020 *
## sp01_env -0.473000 0.0296 -0.4990 -0.3980 *
## sp01_env2 0.276000 0.0687 0.1430 0.4020 *
## sp02_intercept -0.101000 0.0732 -0.2410 0.0477
## sp02_env -0.481000 0.0205 -0.5000 -0.4260 *
## sp02_env2 0.094400 0.0682 -0.0350 0.2330
## sp03_intercept 0.006770 0.0748 -0.1410 0.1520
## sp03_env -0.481000 0.0199 -0.5000 -0.4280 *
## sp03_env2 0.002630 0.0679 -0.1350 0.1300
## sp04_intercept 0.111000 0.0713 -0.0349 0.2420
## sp04_env -0.475000 0.0245 -0.4990 -0.4100 *
## sp04_env2 -0.110000 0.0675 -0.2400 0.0263
## sp05_intercept 0.292000 0.0705 0.1490 0.4300 *
## sp05_env -0.470000 0.0291 -0.4990 -0.3920 *
## sp05_env2 -0.262000 0.0691 -0.3960 -0.1260 *
## sp06_intercept 0.342000 0.0693 0.1990 0.4770 *
## sp06_env -0.463000 0.0339 -0.4990 -0.3750 *
## sp06_env2 -0.329000 0.0655 -0.4550 -0.2050 *
## sp07_intercept 0.469000 0.0285 0.3950 0.4990 *
## sp07_env -0.430000 0.0756 -0.4990 -0.2470 *
## sp07_env2 -0.441000 0.0451 -0.4980 -0.3320 *
## sp08_intercept 0.483000 0.0171 0.4370 0.5000 *
## sp08_env -0.356000 0.0953 -0.4940 -0.1460 *
## sp08_env2 -0.460000 0.0338 -0.4990 -0.3780 *
## sp09_intercept 0.484000 0.0158 0.4390 0.5000 *
## sp09_env -0.362000 0.0874 -0.4920 -0.1740 *
## sp09_env2 -0.478000 0.0211 -0.4990 -0.4220 *
## sp10_intercept 0.488000 0.0117 0.4550 0.5000 *
## sp10_env -0.091000 0.0934 -0.2500 0.1240
## sp10_env2 -0.484000 0.0144 -0.4990 -0.4470 *
## sp11_intercept 0.488000 0.0115 0.4590 0.5000 *
## sp11_env 0.198000 0.0886 0.0421 0.3900 *
## sp11_env2 -0.475000 0.0230 -0.4990 -0.4160 *
## sp12_intercept 0.487000 0.0130 0.4490 0.5000 *
## sp12_env 0.462000 0.0346 0.3710 0.4990 *
## sp12_env2 -0.482000 0.0176 -0.5000 -0.4360 *
## sp13_intercept 0.478000 0.0211 0.4210 0.4990 *
## sp13_env 0.423000 0.0643 0.2730 0.4980 *
## sp13_env2 -0.478000 0.0205 -0.4990 -0.4240 *
## sp14_intercept 0.475000 0.0238 0.4090 0.4990 *
## sp14_env 0.457000 0.0511 0.3140 0.4990 *
## sp14_env2 -0.454000 0.0381 -0.4980 -0.3600 *
## sp15_intercept 0.392000 0.0633 0.2600 0.4920 *
## sp15_env 0.466000 0.0318 0.3810 0.4990 *
## sp15_env2 -0.347000 0.0678 -0.4710 -0.2050 *
## sp16_intercept 0.276000 0.0727 0.1390 0.4220 *
## sp16_env 0.476000 0.0288 0.4080 0.5000 *
## sp16_env2 -0.259000 0.0701 -0.3960 -0.1260 *
## sp17_intercept 0.177000 0.0707 0.0404 0.3200 *
## sp17_env 0.478000 0.0211 0.4210 0.4990 *
## sp17_env2 -0.171000 0.0691 -0.3090 -0.0392 *
## sp18_intercept -0.000261 0.0741 -0.1480 0.1430
## sp18_env 0.481000 0.0197 0.4280 0.4990 *
## sp18_env2 0.009230 0.0702 -0.1350 0.1450
## sp19_intercept -0.092200 0.0739 -0.2300 0.0535
## sp19_env 0.480000 0.0201 0.4240 0.4990 *
## sp19_env2 0.101000 0.0678 -0.0311 0.2390
## sp20_intercept -0.277000 0.0738 -0.4190 -0.1340 *
## sp20_env 0.471000 0.0283 0.3960 0.4990 *
## sp20_env2 0.246000 0.0695 0.1060 0.3820 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 20 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.414, and the DIC is 77410. Computation involved 2500 Gibbs steps, with a burnin of 500. Dimension reduction was implemented with N = 11 and r = 8. Tau <- solve(mod_gjam1$parameters$sigMu)
Tau_n = to_prec(Tau)
#postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
#postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
#pH<-convert_to_m(postH)
#pL<-convert_to_m(postL)
#R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Correlation cor(Y)")
corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("R")
corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("E matrix")
# corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("Tau")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# title("R signif")
corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("True interactions")
data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod)
Environmental filtering + Facilitation dense 5 species
JSDM
mf5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mf5)## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.#mf5$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(mf5$R$B),4)))## Maximum Rhat value for Beta: 1.3081cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(mf5$R$Rho),4)))## Maximum Rhat value for Rho: NAcat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(mf5$R$EnvRho),4)))## Maximum Rhat value for EnvRho: 1.1421cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(mf5$R$Tau),4)))## Warning in max(mf5$R$Tau): no non-missing arguments to max; returning -Inf## Maximum Rhat value for Tau: -Infmf5$n.eff## $B
## [,1] [,2] [,3]
## [1,] 61 47 57
## [2,] 75 496 90
## [3,] 117 305 106
## [4,] 15 15 16
## [5,] 137 161 170
##
## $Rho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 44 44 678 61
## [2,] 44 1 88 199 334
## [3,] 44 88 1 282 110
## [4,] 678 199 282 1 171
## [5,] 61 334 110 171 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 102 59 41 63
## [2,] 102 1 156 33 102
## [3,] 59 156 1 35 225
## [4,] 41 33 35 1 34
## [5,] 63 102 225 34 1mf5$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000data<-sim_data$FacDenseSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm(mf5,data$Y,fac_inter[[4]])
Gjam
data<-sim_data$FacDenseSp5
#fit_gjam(data,5000,500,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
load_gjam(data,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 113.0 61.5 13.9 237
## env2 67.7 31.3 19.6 134
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05
## intercept -1.300 1.020 2.990 1.570 -1.840
## env -1.540 -1.640 0.114 2.570 1.760
## env2 0.558 -0.451 -1.120 0.308 0.299
## RMSPE 0.308 0.315 0.201 0.331 0.308
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -1.300 0.0975 -1.4900 -1.110 *
## sp01_env -1.540 0.1450 -1.8300 -1.300 *
## sp01_env2 0.558 0.0909 0.3960 0.749 *
## sp02_intercept 1.020 0.1180 0.8440 1.300 *
## sp02_env -1.640 0.0966 -1.8200 -1.450 *
## sp02_env2 -0.451 0.0802 -0.6250 -0.298 *
## sp03_intercept 2.990 0.1230 2.7500 3.230 *
## sp03_env 0.114 0.0787 -0.0114 0.287
## sp03_env2 -1.120 0.0715 -1.2600 -0.977 *
## sp04_intercept 1.570 0.0908 1.4100 1.770 *
## sp04_env 2.570 0.1370 2.3300 2.870 *
## sp04_env2 0.308 0.0514 0.2070 0.410 *
## sp05_intercept -1.840 0.1750 -2.1700 -1.510 *
## sp05_env 1.760 0.1520 1.4600 2.060 *
## sp05_env2 0.299 0.0823 0.1450 0.448 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 5 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.296, and the DIC is 80408. Computation involved 5000 Gibbs steps, with a burnin of 500.
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))HMSC
data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, interact = fac_inter[[4]])
Environmental filtering + Facilitation dense 10 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.2866## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.4808## Warning in max(mfd10$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 25 94 22
## [2,] 27 26 25
## [3,] 223 198 394
## [4,] 24 27 23
## [5,] 661 1000 896
## [6,] 22 20 29
## [7,] 275 464 65
## [8,] 445 489 889
## [9,] 16 18 19
## [10,] 19 33 16
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 96 162 47 215 61 204 110 270 93
## [2,] 96 1 200 34 843 55 883 202 223 443
## [3,] 162 200 1 272 1000 801 53 1000 132 757
## [4,] 47 34 272 1 387 126 281 247 83 72
## [5,] 215 843 1000 387 1 233 210 467 590 357
## [6,] 61 55 801 126 233 1 41 33 35 71
## [7,] 204 883 53 281 210 41 1 184 531 55
## [8,] 110 202 1000 247 467 33 184 1 1000 864
## [9,] 270 223 132 83 590 35 531 1000 1 161
## [10,] 93 443 757 72 357 71 55 864 161 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 109 29 42 43 25 42 29 15 15
## [2,] 109 1 62 60 54 33 98 44 12 12
## [3,] 29 62 1 26 580 27 200 1000 17 17
## [4,] 42 60 26 1 30 34 37 32 24 22
## [5,] 43 54 580 30 1 26 163 740 27 34
## [6,] 25 33 27 34 26 1 46 88 26 22
## [7,] 42 98 200 37 163 46 1 499 19 20
## [8,] 29 44 1000 32 740 88 499 1 22 27
## [9,] 15 12 17 24 27 26 19 22 1 105
## [10,] 15 12 17 22 34 22 20 27 105 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
data<-sim_data$FacDenseSp10
#fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
load_gjam(data,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])##
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 160.0 35.60 101.0 240
## env2 18.5 4.88 10.4 29
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08 sp09
## intercept 0.212 -1.490 -0.877 2.410 -0.0951 2.400 0.497 -0.401 1.260
## env -2.760 -2.660 -1.810 -3.010 -0.5910 1.160 1.330 1.030 3.990
## env2 1.070 -0.969 -0.644 0.319 -0.7020 -0.454 -0.256 -0.347 0.095
## RMSPE 0.274 0.375 0.406 0.305 0.4930 0.259 0.390 0.471 0.248
## sp10
## intercept -0.893
## env 2.250
## env2 -0.411
## RMSPE 0.354
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept 0.2120 0.0687 0.0833 0.3520 *
## sp01_env -2.7600 0.1740 -3.1200 -2.4700 *
## sp01_env2 1.0700 0.1350 0.8190 1.3200 *
## sp02_intercept -1.4900 0.2000 -1.9200 -1.1600 *
## sp02_env -2.6600 0.2110 -3.0900 -2.3000 *
## sp02_env2 -0.9690 0.1130 -1.2400 -0.7780 *
## sp03_intercept -0.8770 0.1070 -1.0900 -0.6900 *
## sp03_env -1.8100 0.1470 -2.1200 -1.5700 *
## sp03_env2 -0.6440 0.1190 -0.9250 -0.4570 *
## sp04_intercept 2.4100 0.2040 2.0500 2.8000 *
## sp04_env -3.0100 0.2470 -3.4600 -2.5900 *
## sp04_env2 0.3190 0.0696 0.1860 0.4540 *
## sp05_intercept -0.0951 0.0554 -0.2050 0.0147
## sp05_env -0.5910 0.0628 -0.7140 -0.4700 *
## sp05_env2 -0.7020 0.0743 -0.8490 -0.5590 *
## sp06_intercept 2.4000 0.1260 2.1600 2.6400 *
## sp06_env 1.1600 0.0761 1.0100 1.3100 *
## sp06_env2 -0.4540 0.0649 -0.5840 -0.3310 *
## sp07_intercept 0.4970 0.0563 0.3840 0.6040 *
## sp07_env 1.3300 0.0797 1.1800 1.4900 *
## sp07_env2 -0.2560 0.0554 -0.3670 -0.1470 *
## sp08_intercept -0.4010 0.0539 -0.5070 -0.2950 *
## sp08_env 1.0300 0.0877 0.8540 1.1900 *
## sp08_env2 -0.3470 0.0592 -0.4640 -0.2320 *
## sp09_intercept 1.2600 0.0922 1.0900 1.4500 *
## sp09_env 3.9900 0.3810 3.3800 4.6800 *
## sp09_env2 0.0950 0.0613 -0.0361 0.2050
## sp10_intercept -0.8930 0.0796 -1.0700 -0.7550 *
## sp10_env 2.2500 0.1280 2.0400 2.5200 *
## sp10_env2 -0.4110 0.0576 -0.5230 -0.2980 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 10 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.367, and the DIC is 157914. Computation involved 5000 Gibbs steps, with a burnin of 500.
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))HMSC
Environmental filtering + Facilitation dense 20 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.6132## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.508## Warning in max(mfd20$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 13 57 12
## [2,] 46 42 42
## [3,] 301 176 90
## [4,] 10 29 10
## [5,] 46 37 35
## [6,] 34 32 37
## [7,] 536 181 102
## [8,] 868 1000 388
## [9,] 121 172 238
## [10,] 100 625 303
## [11,] 499 1000 121
## [12,] 1000 184 395
## [13,] 269 166 153
## [14,] 54 149 179
## [15,] 112 231 138
## [16,] 108 227 63
## [17,] 42 117 37
## [18,] 44 43 71
## [19,] 16 29 14
## [20,] 23 31 27
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 353 436 46 381 42 169 83 93 24 171 130 160
## [2,] 353 1 392 339 33 80 150 134 41 26 184 175 106
## [3,] 436 392 1 51 59 101 171 59 297 569 166 359 111
## [4,] 46 339 51 1 1000 90 35 68 27 118 78 50 40
## [5,] 381 33 59 1000 1 269 156 867 283 48 130 61 46
## [6,] 42 80 101 90 269 1 33 79 162 336 44 51 220
## [7,] 169 150 171 35 156 33 1 303 143 127 97 190 66
## [8,] 83 134 59 68 867 79 303 1 147 109 472 344 52
## [9,] 93 41 297 27 283 162 143 147 1 292 314 137 62
## [10,] 24 26 569 118 48 336 127 109 292 1 367 362 180
## [11,] 171 184 166 78 130 44 97 472 314 367 1 326 229
## [12,] 130 175 359 50 61 51 190 344 137 362 326 1 1000
## [13,] 160 106 111 40 46 220 66 52 62 180 229 1000 1
## [14,] 243 134 169 37 160 166 73 44 36 131 100 182 173
## [15,] 113 82 622 46 37 51 136 145 72 46 273 159 85
## [16,] 555 52 111 58 47 64 285 64 22 28 151 525 36
## [17,] 360 84 62 79 73 569 191 106 45 52 120 102 29
## [18,] 129 128 86 191 621 197 226 202 257 53 32 1000 817
## [19,] 119 73 102 119 348 70 168 41 35 163 41 52 29
## [20,] 190 194 199 117 127 376 84 122 674 46 124 70 40
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 243 113 555 360 129 119 190
## [2,] 134 82 52 84 128 73 194
## [3,] 169 622 111 62 86 102 199
## [4,] 37 46 58 79 191 119 117
## [5,] 160 37 47 73 621 348 127
## [6,] 166 51 64 569 197 70 376
## [7,] 73 136 285 191 226 168 84
## [8,] 44 145 64 106 202 41 122
## [9,] 36 72 22 45 257 35 674
## [10,] 131 46 28 52 53 163 46
## [11,] 100 273 151 120 32 41 124
## [12,] 182 159 525 102 1000 52 70
## [13,] 173 85 36 29 817 29 40
## [14,] 1 90 96 127 116 154 43
## [15,] 90 1 127 137 82 28 237
## [16,] 96 127 1 279 122 95 192
## [17,] 127 137 279 1 136 151 227
## [18,] 116 82 122 136 1 90 236
## [19,] 154 28 95 151 90 1 57
## [20,] 43 237 192 227 236 57 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 16 18 17 18 15 14 15 15 17 18 17 15
## [2,] 16 1 105 12 62 36 45 71 75 190 103 139 74
## [3,] 18 105 1 17 87 57 112 130 220 93 266 109 139
## [4,] 17 12 17 1 11 17 12 12 12 12 15 13 12
## [5,] 18 62 87 11 1 43 62 92 135 162 193 102 124
## [6,] 15 36 57 17 43 1 48 39 39 32 47 33 35
## [7,] 14 45 112 12 62 48 1 551 219 209 874 259 208
## [8,] 15 71 130 12 92 39 551 1 755 425 1000 409 416
## [9,] 15 75 220 12 135 39 219 755 1 407 1000 501 1000
## [10,] 17 190 93 12 162 32 209 425 407 1 668 1000 340
## [11,] 18 103 266 15 193 47 874 1000 1000 668 1 708 1000
## [12,] 17 139 109 13 102 33 259 409 501 1000 708 1 554
## [13,] 15 74 139 12 124 35 208 416 1000 340 1000 554 1
## [14,] 15 78 76 12 133 36 314 516 264 557 591 322 306
## [15,] 14 72 147 12 83 32 133 231 1000 222 611 434 1000
## [16,] 15 26 138 15 40 93 84 70 78 53 119 59 78
## [17,] 80 89 71 23 76 41 35 42 64 55 57 57 50
## [18,] 19 44 92 16 39 44 230 151 136 127 223 241 132
## [19,] 32 22 17 38 21 23 17 17 18 22 22 23 18
## [20,] 21 78 48 17 40 21 27 31 34 39 38 34 30
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 15 14 15 80 19 32 21
## [2,] 78 72 26 89 44 22 78
## [3,] 76 147 138 71 92 17 48
## [4,] 12 12 15 23 16 38 17
## [5,] 133 83 40 76 39 21 40
## [6,] 36 32 93 41 44 23 21
## [7,] 314 133 84 35 230 17 27
## [8,] 516 231 70 42 151 17 31
## [9,] 264 1000 78 64 136 18 34
## [10,] 557 222 53 55 127 22 39
## [11,] 591 611 119 57 223 22 38
## [12,] 322 434 59 57 241 23 34
## [13,] 306 1000 78 50 132 18 30
## [14,] 1 248 61 40 106 21 29
## [15,] 248 1 63 54 124 18 30
## [16,] 61 63 1 29 124 14 24
## [17,] 40 54 29 1 49 59 112
## [18,] 106 124 124 49 1 23 25
## [19,] 21 18 14 59 23 1 42
## [20,] 29 30 24 112 25 42 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Gjam
data<-sim_data$FacDenseSp20
fit_gjam(data,5000,500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])##
## Observations and responses:
## [1] 500 20
## ===========================================================================## expanding covariance chains## ===========================================================================
## Sensitivity by predictor variables f:
## Estimate SE CI_025 CI_975
## env 307 54.9 217 427
## env2 187 34.7 122 262
##
## Coefficient matrix B:
## sp01 sp02 sp03 sp04 sp05 sp06 sp07 sp08
## intercept -0.480 -0.99800 -0.423 -0.180 -0.0128 0.385 0.0559 0.244
## env -1.980 -1.74000 -1.860 -2.680 -1.3400 -1.960 -0.9680 -0.722
## env2 0.308 0.00803 -0.255 -0.555 -0.4870 -0.893 -0.6740 -0.899
## RMSPE 0.321 0.35800 0.370 0.296 0.4070 0.291 0.4510 0.445
## sp09 sp10 sp11 sp12 sp13 sp14 sp15 sp16 sp17
## intercept 0.640 0.712 0.317 0.262 0.579 0.250 0.0141 -0.209 -0.0225
## env -0.914 -0.402 0.105 0.469 0.932 1.050 1.5200 2.100 2.1400
## env2 -1.250 -1.120 -0.865 -0.937 -0.874 -0.957 -1.0500 -0.612 -0.1930
## RMSPE 0.383 0.413 0.471 0.453 0.395 0.407 0.3710 0.331 0.3260
## sp18 sp19 sp20
## intercept 0.0822 0.125 -0.385
## env 2.4300 2.800 2.080
## env2 0.2020 0.794 0.260
## RMSPE 0.2980 0.278 0.317
##
## Coefficient matrix B:
## Estimate SE CI_025 CI_975 sig95
## sp01_intercept -0.48000 0.0709 -0.61400 -0.3370 *
## sp01_env -1.98000 0.1840 -2.28000 -1.6500 *
## sp01_env2 0.30800 0.0638 0.18000 0.4320 *
## sp02_intercept -0.99800 0.1380 -1.27000 -0.7540 *
## sp02_env -1.74000 0.1240 -1.99000 -1.5200 *
## sp02_env2 0.00803 0.0846 -0.16100 0.1510
## sp03_intercept -0.42300 0.1680 -0.75800 -0.2050 *
## sp03_env -1.86000 0.1500 -2.12000 -1.5500 *
## sp03_env2 -0.25500 0.1150 -0.45600 -0.0566 *
## sp04_intercept -0.18000 0.0555 -0.28400 -0.0712 *
## sp04_env -2.68000 0.1380 -2.95000 -2.4200 *
## sp04_env2 -0.55500 0.0641 -0.67700 -0.4330 *
## sp05_intercept -0.01280 0.0903 -0.17800 0.1530
## sp05_env -1.34000 0.0753 -1.49000 -1.2000 *
## sp05_env2 -0.48700 0.1010 -0.65700 -0.2960 *
## sp06_intercept 0.38500 0.0585 0.27000 0.4990 *
## sp06_env -1.96000 0.1200 -2.21000 -1.7400 *
## sp06_env2 -0.89300 0.0928 -1.06000 -0.6990 *
## sp07_intercept 0.05590 0.0662 -0.07310 0.1860
## sp07_env -0.96800 0.0724 -1.12000 -0.8310 *
## sp07_env2 -0.67400 0.0632 -0.80400 -0.5520 *
## sp08_intercept 0.24400 0.0603 0.13200 0.3630 *
## sp08_env -0.72200 0.0649 -0.84700 -0.5930 *
## sp08_env2 -0.89900 0.0723 -1.05000 -0.7690 *
## sp09_intercept 0.64000 0.0555 0.53100 0.7490 *
## sp09_env -0.91400 0.0596 -1.03000 -0.7970 *
## sp09_env2 -1.25000 0.0745 -1.40000 -1.1100 *
## sp10_intercept 0.71200 0.0581 0.59800 0.8230 *
## sp10_env -0.40200 0.0555 -0.51100 -0.2930 *
## sp10_env2 -1.12000 0.0838 -1.27000 -0.9570 *
## sp11_intercept 0.31700 0.0807 0.13400 0.4750 *
## sp11_env 0.10500 0.0607 -0.01620 0.2190
## sp11_env2 -0.86500 0.0678 -0.99500 -0.7330 *
## sp12_intercept 0.26200 0.0541 0.15600 0.3680 *
## sp12_env 0.46900 0.0547 0.36200 0.5740 *
## sp12_env2 -0.93700 0.0660 -1.06000 -0.8050 *
## sp13_intercept 0.57900 0.0808 0.41100 0.7210 *
## sp13_env 0.93200 0.0628 0.81000 1.0500 *
## sp13_env2 -0.87400 0.0659 -1.01000 -0.7510 *
## sp14_intercept 0.25000 0.0625 0.13100 0.3720 *
## sp14_env 1.05000 0.0695 0.90700 1.1800 *
## sp14_env2 -0.95700 0.1110 -1.14000 -0.7350 *
## sp15_intercept 0.01410 0.0589 -0.09820 0.1310
## sp15_env 1.52000 0.0761 1.38000 1.6700 *
## sp15_env2 -1.05000 0.0809 -1.21000 -0.9040 *
## sp16_intercept -0.20900 0.0883 -0.37100 -0.0447 *
## sp16_env 2.10000 0.1020 1.91000 2.3100 *
## sp16_env2 -0.61200 0.0559 -0.72300 -0.5050 *
## sp17_intercept -0.02250 0.0508 -0.12200 0.0783
## sp17_env 2.14000 0.1300 1.92000 2.4000 *
## sp17_env2 -0.19300 0.0816 -0.34100 -0.0369 *
## sp18_intercept 0.08220 0.0547 -0.02540 0.1870
## sp18_env 2.43000 0.1090 2.22000 2.6500 *
## sp18_env2 0.20200 0.0628 0.07750 0.3160 *
## sp19_intercept 0.12500 0.0624 0.00441 0.2470 *
## sp19_env 2.80000 0.1190 2.57000 3.0300 *
## sp19_env2 0.79400 0.0876 0.63200 0.9540 *
## sp20_intercept -0.38500 0.0727 -0.51200 -0.2370 *
## sp20_env 2.08000 0.0929 1.90000 2.2600 *
## sp20_env2 0.26000 0.0780 0.10300 0.3960 *
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Coefficient matrix B, standardized for X and W:
## NULL
##
## Last column indicates if 95% posterior distribution contains zero.
##
## Design Table
## env env2
## VIF 1 1
## env2 0 NA
##
## Sample contains n = 500 observations on S = 20 response variables, and 2 predictors. Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.374, and the DIC is 395024. Computation involved 5000 Gibbs steps, with a burnin of 500.
#load_gjam(data,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))HMSC
data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, interact = fac_inter[[6]])
Environmental filtering + Facilitation sparse 5 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.373## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.4048## Warning in max(mfs5$R$Tau): no non-missing arguments to max; returning -Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 17 44 17
## [2,] 169 469 127
## [3,] 1000 1000 500
## [4,] 20 25 20
## [5,] 13 66 14
##
## $Rho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 111 92 123 819
## [2,] 111 1 283 40 143
## [3,] 92 283 1 205 28
## [4,] 123 40 205 1 102
## [5,] 819 143 28 102 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 20 19 29 13
## [2,] 20 1 252 20 18
## [3,] 19 252 1 21 17
## [4,] 29 20 21 1 32
## [5,] 13 18 17 32 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Environmental filtering + Facilitation sparse 10 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.5479## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.4746## Warning in max(mfs10$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 52 49 55
## [2,] 11 17 10
## [3,] 19 20 24
## [4,] 1000 374 266
## [5,] 1000 567 333
## [6,] 262 286 296
## [7,] 839 602 1000
## [8,] 129 212 94
## [9,] 462 346 407
## [10,] 49 251 47
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 122 94 384 296 316 114 140 292 337
## [2,] 122 1 558 143 481 155 81 262 144 148
## [3,] 94 558 1 89 355 800 1000 317 214 104
## [4,] 384 143 89 1 172 290 191 213 543 438
## [5,] 296 481 355 172 1 273 110 421 192 116
## [6,] 316 155 800 290 273 1 621 178 401 101
## [7,] 114 81 1000 191 110 621 1 179 140 53
## [8,] 140 262 317 213 421 178 179 1 630 181
## [9,] 292 144 214 543 192 401 140 630 1 213
## [10,] 337 148 104 438 116 101 53 181 213 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 15 31 147 111 128 138 80 119 41
## [2,] 15 1 37 13 14 15 12 22 11 15
## [3,] 31 37 1 26 25 30 25 30 23 97
## [4,] 147 13 26 1 449 363 1000 78 1000 48
## [5,] 111 14 25 449 1 615 361 97 630 70
## [6,] 128 15 30 363 615 1 696 109 358 50
## [7,] 138 12 25 1000 361 696 1 110 488 48
## [8,] 80 22 30 78 97 109 110 1 153 34
## [9,] 119 11 23 1000 630 358 488 153 1 68
## [10,] 41 15 97 48 70 50 48 34 68 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Environmental filtering + Facilitation sparse 20 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.5653## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.4826## Warning in max(mfs20$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 24 21 32
## [2,] 22 59 19
## [3,] 27 108 32
## [4,] 22 251 24
## [5,] 227 23 216
## [6,] 273 395 140
## [7,] 351 1000 250
## [8,] 177 104 232
## [9,] 206 144 69
## [10,] 384 290 280
## [11,] 128 404 1000
## [12,] 1000 570 505
## [13,] 887 747 270
## [14,] 170 1000 1000
## [15,] 175 51 217
## [16,] 43 147 29
## [17,] 10 29 13
## [18,] 22 26 16
## [19,] 12 84 13
## [20,] 23 33 25
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 422 57 96 17 43 24 74 81 32 502 101 59
## [2,] 422 1 91 212 54 128 110 167 93 75 89 128 87
## [3,] 57 91 1 172 94 434 162 154 111 298 118 29 33
## [4,] 96 212 172 1 92 399 138 52 123 41 74 65 43
## [5,] 17 54 94 92 1 164 95 612 312 19 71 18 95
## [6,] 43 128 434 399 164 1 152 162 180 349 118 37 80
## [7,] 24 110 162 138 95 152 1 379 185 129 435 48 101
## [8,] 74 167 154 52 612 162 379 1 105 519 122 477 204
## [9,] 81 93 111 123 312 180 185 105 1 112 61 263 190
## [10,] 32 75 298 41 19 349 129 519 112 1 247 184 102
## [11,] 502 89 118 74 71 118 435 122 61 247 1 120 301
## [12,] 101 128 29 65 18 37 48 477 263 184 120 1 92
## [13,] 59 87 33 43 95 80 101 204 190 102 301 92 1
## [14,] 150 38 72 603 57 373 83 739 91 354 80 100 193
## [15,] 45 80 61 273 90 38 24 82 39 153 99 95 137
## [16,] 92 111 36 23 43 83 107 95 198 34 55 302 66
## [17,] 97 28 53 253 137 244 115 38 170 26 25 104 118
## [18,] 107 38 109 107 103 72 57 73 98 121 44 47 27
## [19,] 46 150 110 517 52 216 66 81 42 23 137 80 54
## [20,] 51 66 339 154 64 75 85 98 68 98 82 40 31
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 150 45 92 97 107 46 51
## [2,] 38 80 111 28 38 150 66
## [3,] 72 61 36 53 109 110 339
## [4,] 603 273 23 253 107 517 154
## [5,] 57 90 43 137 103 52 64
## [6,] 373 38 83 244 72 216 75
## [7,] 83 24 107 115 57 66 85
## [8,] 739 82 95 38 73 81 98
## [9,] 91 39 198 170 98 42 68
## [10,] 354 153 34 26 121 23 98
## [11,] 80 99 55 25 44 137 82
## [12,] 100 95 302 104 47 80 40
## [13,] 193 137 66 118 27 54 31
## [14,] 1 69 337 73 53 44 61
## [15,] 69 1 88 19 33 36 34
## [16,] 337 88 1 210 83 48 51
## [17,] 73 19 210 1 51 125 51
## [18,] 53 33 83 51 1 123 37
## [19,] 44 36 48 125 123 1 37
## [20,] 61 34 51 51 37 37 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## [1,] 1 42 28 93 37 48 50 58 36 53 87 66 47
## [2,] 42 1 15 128 56 25 37 24 36 35 29 28 29
## [3,] 28 15 1 29 57 32 42 62 37 51 84 56 50
## [4,] 93 128 29 1 25 32 26 32 37 50 39 31 33
## [5,] 37 56 57 25 1 137 521 146 787 251 123 139 141
## [6,] 48 25 32 32 137 1 536 243 205 335 240 303 169
## [7,] 50 37 42 26 521 536 1 189 402 349 214 254 319
## [8,] 58 24 62 32 146 243 189 1 251 258 1000 1000 283
## [9,] 36 36 37 37 787 205 402 251 1 705 162 189 266
## [10,] 53 35 51 50 251 335 349 258 705 1 244 246 1000
## [11,] 87 29 84 39 123 240 214 1000 162 244 1 1000 367
## [12,] 66 28 56 31 139 303 254 1000 189 246 1000 1 324
## [13,] 47 29 50 33 141 169 319 283 266 1000 367 324 1
## [14,] 36 28 50 43 160 187 217 312 412 1000 276 246 749
## [15,] 55 31 48 36 89 301 132 550 140 294 1000 1000 305
## [16,] 72 22 40 57 26 43 37 50 46 76 73 56 58
## [17,] 17 47 24 37 15 15 14 20 15 17 26 23 15
## [18,] 76 74 21 51 23 30 21 33 22 29 51 36 23
## [19,] 108 155 40 28 32 17 19 19 24 22 21 18 18
## [20,] 41 182 143 38 68 42 39 31 55 57 34 33 38
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 36 55 72 17 76 108 41
## [2,] 28 31 22 47 74 155 182
## [3,] 50 48 40 24 21 40 143
## [4,] 43 36 57 37 51 28 38
## [5,] 160 89 26 15 23 32 68
## [6,] 187 301 43 15 30 17 42
## [7,] 217 132 37 14 21 19 39
## [8,] 312 550 50 20 33 19 31
## [9,] 412 140 46 15 22 24 55
## [10,] 1000 294 76 17 29 22 57
## [11,] 276 1000 73 26 51 21 34
## [12,] 246 1000 56 23 36 18 33
## [13,] 749 305 58 15 23 18 38
## [14,] 1 255 75 16 24 19 39
## [15,] 255 1 125 19 38 15 36
## [16,] 75 125 1 13 48 14 30
## [17,] 16 19 13 1 103 14 14
## [18,] 24 38 48 103 1 18 25
## [19,] 19 15 14 14 18 1 21
## [20,] 39 36 30 14 25 21 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Environmental filtering + Competition dense 5 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## Successful convergence based on Rhat values (all < 1.1).## Maximum Rhat value for Beta: 1.0324## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.0313## Warning in max(mcmpd5$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 1000 1000 389
## [2,] 210 355 355
## [3,] 320 1000 272
## [4,] 1000 1000 1000
## [5,] 1000 463 95
##
## $Rho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 252 209 114 59
## [2,] 252 1 132 483 52
## [3,] 209 132 1 611 96
## [4,] 114 483 611 1 162
## [5,] 59 52 96 162 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 556 455 431 207
## [2,] 556 1 795 661 120
## [3,] 455 795 1 1000 130
## [4,] 431 661 1000 1 100
## [5,] 207 120 130 100 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Environmental filtering + Competition dense 10 species
## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.1845## Maximum Rhat value for Rho: NA## Maximum Rhat value for EnvRho: 1.1771## Warning in max(mcmpd10$R$Tau): no non-missing arguments to max; returning -
## Inf## Maximum Rhat value for Tau: -Inf## $B
## [,1] [,2] [,3]
## [1,] 24 38 25
## [2,] 802 157 93
## [3,] 833 870 1000
## [4,] 573 462 566
## [5,] 994 601 648
## [6,] 708 1000 329
## [7,] 420 240 320
## [8,] 29 28 33
## [9,] 123 174 114
## [10,] 206 117 100
##
## $Rho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 49 197 83 87 125 89 176 184 124
## [2,] 49 1 40 82 315 39 73 187 71 174
## [3,] 197 40 1 387 377 273 177 209 857 269
## [4,] 83 82 387 1 81 136 1000 50 131 442
## [5,] 87 315 377 81 1 1000 118 169 59 219
## [6,] 125 39 273 136 1000 1 1000 144 242 144
## [7,] 89 73 177 1000 118 1000 1 55 216 41
## [8,] 176 187 209 50 169 144 55 1 381 155
## [9,] 184 71 857 131 59 242 216 381 1 175
## [10,] 124 174 269 442 219 144 41 155 175 1
##
## $EnvRho
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 39 33 59 59 34 24 134 31 59
## [2,] 39 1 96 83 84 108 170 178 92 306
## [3,] 33 96 1 509 591 612 105 37 105 96
## [4,] 59 83 509 1 1000 326 105 44 114 97
## [5,] 59 84 591 1000 1 385 119 39 181 99
## [6,] 34 108 612 326 385 1 256 59 112 172
## [7,] 24 170 105 105 119 256 1 113 67 194
## [8,] 134 178 37 44 39 59 113 1 69 222
## [9,] 31 92 105 114 181 112 67 69 1 123
## [10,] 59 306 96 97 99 172 194 222 123 1## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000
Environmental filtering + Competition dense 20 species
Environmental filtering + Competition sparse 5 species
Environmental filtering + Competition dense 10 species
Environmental filtering + Competition dense 5 species
Graph
Visualization for true interactions
